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RESUMO 


Neste trabalho, após um repasse dos distintos métodos de 
cálculo numérico dos factores de intensificacao de tensoes, apresenta se 
um novo elemento finito singular que calcula -se baíxo condicoes de carga 
mistas (modo I e II combinados). Como conclusao extrain-se algumas con- 
sequéncias sobre o uso deste elemento e dans-se algumas recomendacoes para 
a sua utilizacao. 


ABSTRACT 


In the present paper after reviewíng the different methods for 
evaluating the stress intensity factors, we present a new singular finite 
element, which is evaluated under missed loading conditions (Mode I and 
Mode II are presented). As conclusion there are taken out some results 
about the performance of these element, and some advices for its use are 
given. 


RESUME 


Dans ce travail, apres faire una revision des differentes 
methodes de calcul numérique des facteurs d'intensité des contraintes, on 
présente un nouvel elément finí, quí est evalué sous conditions míxtes de 
chargement (mode I et mode Il combinés). On concluit quelques conséquences 
sur le performance de cet element, est on fait quelques recommendations 
pour son emploi. 
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1. INTRODUCCION 


Como es bien sabido las soluciones analíticas de problemas de 
mecánica de la fractura se limitan a situaciones idealizadas donde: (i) el 
dominio se considera infinito, (ii) el material es homogéneo e isótropo y 
(iii) las condiciones de contorno son relativamente simples. Para poder 
abordar problemas prácticos sobre estructuras de cualquier forma, con- 
diciones de contorno complejas y propiedades arbitrarias, se hacen im- 
prescindibles las técnicas numéricas. En este sentido el método de los 
elementos finitos es una herramienta de uso general y muy extendido en el 
campo de la mecénica de la fractura. 


Existen varios elementos finitos especiales para problemas de 
fractura plana [1] que intentan reproducir el comportamiento singular de 
la solución (tensiones o deformaciones) en el vértice de la fisura. Todos 
ellos han sido abundantemente analizados y comprobados en la bibliografía 
(ver por ejemplo [1], [2]), pero en general se observaba una laguna en el 
análisis teórico de la solución aproximada dentro del elemento. Este 
estudio ha sido abordado recientemente [3], [4] y muestra claramente las 
razones por las cuales algunos elementos finitos producían resultados 
anómalos en algunos casos. Por otro lado se observa claramente que la 
mayoría de los tests numéricos efectuados han sido llevados a cabo con- 
siderando condiciones de carga simples (concretamente en modo 1). Sin em- 
bargo es sabido que la prueba más severa para elementos singulares es en 
modo mixto donde están presentes el Modo 1 y el Modo 11. En el presente 
trabajo se pretende efectuar una evaluación numérica de un nuevo elemento 
finito singular, recientemente introducido [5], [6], bajo condiciones de 
carga combinadas. Para ello primeramente se harán algunas consideraciones 
sobre el cálculo numérico de los factores de intensificación de tensiones, 
que en cierto modo justificarán la solución propuesta (82). En el $3 se 
formula un nuevo elemento finito singular isoparamétrico de lados curvos, 
y en el 84 se considera una pieza sometida a una carga en modo mixto, que 
se analiza con el método propuesto haciéndose varias recomendaciones sobre 
el uso de elementos singulares. 


2. “SOBRE EL CALCULO NUMERICO DE LOS FACTORES DE INTENSIFICACION DE 
TENSIONES 


Consideremos un sólido fisurado sometido a la acción de una 
serie de fuerzas exteriores. La resolución de las ecuaciones de equilibrio 
cuasi-estático en un entorno del vértice de la grieta y bajo las hipótesis 
de que aquél se encuentra en un medio linealmente elástico, isótropo y 
homogéneo conduce a las siguientes componentes asintóticas del campo de 
desplazamientos y tensiones: 


1 1 
ui E Ki ñ 2u Jr/2x Ei (9) + Kr1 2u J1/27 8i (6) i=1,2 
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. 7 83 (0 i=1,2 1b 
77 y O 1 (1b) 


1 
043 = K R=m— ii (0) +K 
ij 1 (nr) ij 11 
En las expresiones anteriores los subíndices hacen referencia 
al sistema de ejes de la figura 1, r y 9 son las coordenadas polares 
centradas en S (vértice de la grieta). 


Las funciones angulares f,, B;, f., Y Bug sólo dependen del 
ángulo $ y de la aproximación plana elegida, p es el módulo de cizal- 
lamiento y K, y K,, los factores de intensificación de tensiones en modo 1 
y II. 


La expresión (1b) indica que las tensiones se hacen infinitas 
en el vértice dela grieta. Físicamente existe sin embargo una relajación 
de las tensiones debido a deformaciones inelásticas y plastificación en un 
entorno de S. Si la zona donde se producen estos fenómenos es pequeña com- 
parada con la singular (donde tienen aplicación las (1)) podemos tratar el 
problema como elástico, lineal, (mecánica de la fractura lineal y 
elástica, MFLE). Restringiéndonos por tanto a la MFLE podemos decir que 
todos los criterios de fractura se reducen a una relación del tipo: 

Í(K 


K constantes del material) = O (2) 


E 11? 


Es por tanto fundamental según indica (2) el cálculo de los 
factores K, y K,,. Para ello se puede comenzar analizando las expresiones 
(1). Si la solución del problema elástico es conocida (mediante un 
análisis por el método de los elementos finitos por ejemplo), Kiss Y KR; 
pueden ser evaluados en cualquier punto próximo a S donde son válidas las 
expresiones asintóticas (1). Haciendo dicho cálculo para 0 = rm y 0 = -1 se 
llega a las siguientes expresiones: 


Ki = lim J(2r1) 029 (x1 PS r,0) Es 


r>0 

2u 
= — lim /2n/r (uy (x1=r7 ,+0) - us (Xx =r, -0)) 

k+1 r>0 


o lim /(2rr) 012 (1 = r,0) = 


r>0 
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2p 
== im 2/0 (01 59,470) 1 (0-0)) (3a,b) 
k+1 1>0 


Las expresiones (3) indican que los factores de in- 
tensificación pueden interpretarse como factores de discontinuidad de los 
desplazamientos [7] y [8], y en este sentido, si mediante un cálculo por 
el método de elementos finitos conocemos el campo de desplazamientos en 
puntos "suficientemente" próximos a S (r>0), se llega a las expresiones 
(figura 2). 


2p 
K+ = — /2n/t (uo% - usd) 
1 1 2 2 


2u A B 
K11 = ce ¿2n/x Cu" 11 (4a,b) 


donde u,*,u,*,u,*,u,? son los desplazamientos verticales y horizontales de 


los nodos A y B (inicialmente superpuestos). 


Las expresiones (4) constituyen la llamada técnica de cor- 
relación de los desplazamientos y dan buenos resultados siempre que A y B 
se encuentren suficientemente cerca del vértice de grieta, [9]. 


A partir de las soluciones análíticas (1), se pueden obtener 
nuevas aproximaciones de los factores de intensificación de tensiones 
cuando la zona singular es modelizada mediante elementos finitos 
singulares o especiales. De entre éstos los más populares son los 
denominados "a 1/4", para los cuales se tiene un valor mejorado de K, y K,, 
mediante la fórmula: 


K*, = 2 K,, - L 


A 4v i=IóÓ II (5) 


donde K, es el factor evaluado en el "medio" del elemento (r = 1/4) y K, 
en el vértice de dicho elemento (r=1). Las (5) desarrolladas conducen a 
las siguientes expresiones para K, y K,,, [10]: 





di (l+1) Inf21 . (u3%-u95) + (u3P-uz) 
2 AB D, E 


(k+1) 


Expresiones análogas a la (5) para elementos singulares de 
tipo Akin [11], han sido obtenidos en modo 1 [10], en concreto un valor 
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mejorado de K, es: 


A (7) 


En definitiva se observa claramente que mediante un cálculo 
por el método de elementos finitos pueden obtenerse en una etapa de post- 
proceso (pero sin ningún esfuerzo adicional) los factores K, y K,, ya sea 
por (4), (5), (6) ó (7). La exactitud de los resultados sólo depende de la 
capacidad de los elementos usados para aproximar los campos asintóticos 


(dy. 


Existen, sin embargo, otras técnicas bien conocidas para la 
evaluación de K, y K,,, que podríamos llamar energéticas. Una de ellas 
sería calcular, mediante una integral de contorno la integral J de Rice, 
(que en el marco de la MFLE coincide con la tasa de restitución de ener- 
gía). Esta integral permite calcular K,? 6 K,,? cuando la carga y la 
geometría son tales que se pueden encontrar los modos puramente simétricos 
o puramente antisimétricos. Sin embargo en modo mixto J nos da la forma 
cuadrática K,? + As Una forma de resolver el problema del desacopla- 
miento entre el modo 1 y 1I es calcular las integrales asociadas J, y J;y,, 
[12], ello exige calcular las partes pares e impares de las componentes de 
u respecto a Xx, y posteriormente dos integrales de contorno extendidas 
sobre un dominio simétrico respecto a Ox,. 


Recientemente se ha propuesto un nuevo método de cálculo de la 
tasa de restitución de energía, G, [13]. Este método, que ha ya sido im- 
plementado en un gran código de cálculo, [14], permite calcular G mediante 
la expresión: 


ón a ón a 
A A 


(8) 
OX t dy t 


En la expresión (8) rm es el llamado indicador variacional, que 
es la energía potencial total del sólido fracturado, t, es el espesor de 
la pieza y 91/0x, ó0r/0y son las integrales asociadas J; y J¡, de donde se 
obtienen K, y K,,. Desgraciadamente (8) exige el conocimiento de a,, a,, 
que son las componentes de un vector unitario en la dirección en que se 
propagaría la grieta. Esta dirección es sin embargo desconocida a priori, 
si se consideran como válidos los criterios de fractura en modo mixto S y 
T, comúnmente admitidos, que exigen, para su formulación el conocimiento 
de E; y Ky 115]. 


Es claro por tanto que ambos métodos energéticos exigen una 
etapa de postproceso elaborada, y son por tanto computacionalmente más 
costosos que los métodos basados en (1). Para estos últimos se requieren 
sin embargo elementos finitos capaces de reproducir el comportamiento 
asintótico de la solución en S. Un tal elemento se propone a continuación. 
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3. UN_ELEMENTO FINITO SINGULAR 


Desde hace tiempo [16] es sabido que una forma simple de 
generar elementos finitos singulares en MFLE consiste en usar elementos 
cuadráticos triangulares (normales o degenerados) o cuadrilaterales y des- 
plazar los nodos intermedios más próximos al vértice de grieta a una dis- 
tancia igual a 1/4 de la longitud del lado correspondiente. 


Esta técnica tan simple ha dado resultados muy satisfactorios 
en problemas simples (principalmente en modo 1). Posteriormente un 
análisis más detallado de dichos elementos ha sido llevado a cabo, al ob- 
jeto de estudiar el comportamiento de la solución aproximada a lo largo de 
líneas radiales [8], [17]. Para lo que nos ocupa las conclusiones más im- 
portantes son: 


i) los cuadriláteros no presentan un buen comportamiento, salvo en los 
propios lados del elemento, 


ii) todos los elementos deben tener los lados rectos. 


Se pueden relajar estos dos restricciones, es decir constuir 
elementos cuadrilaterales de lados curvos adecuados para la MFLE, 
basándose en la técnica propuesta po J. Aalto [5]. 


Como es bien sabido cuando se trabaja con elementos 
isoparamétricos cuadráticos todos los cálculos a nivel elemental son 
llevados a cabo sobre un elemento de referencia ubicado en un plano 
paramétrico de coordenadas que denominares (€, n), -1< €, n < 1. Sea T la 
transformación geométrica que aplica dicho elemento sobre el elemento real 
ubicado en el plano (x,y). Por otro lado el método de los elementos 
finitos consiste en construir subespacios de dimensión finita de espacios 
vectoriales de dimensión infinita mediante la técnica de interpolación 
polinomial a trozos, tomando como soporte de interpelación los nodos de 
los elementos; sea rn la interpolación elegida, es decir la aplicación que 
para cada punto (€, y) nos dá el valor aproximado de la solución. Por otro 
lado alrededor de los puntos singulares es conocida la solución analítica, 
mediante un desarrollo asintótico (1). Sea A dicha solución. 


Por motivos de simplicidad vamos a suponer que el tema que nos 
ocupa es un problema de flujo de agua a través de un medio poroso, donde 
las funciones a aproximar son el potencial de velocidad $ y el potencial 
vector de velocidad v (suponemos movimiento plano). Tendremos entonces el 
siguiente diagrama de aplicaciones: 


T A 
(€, 1) > (2,y) > (0,0) (9) 
| 





T t 


Es decir, 7 = A o T 
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Con relación a la figura 3, consideremos el punto singular 
S(2r-a>r) y una zona próxima a S limitada por líneas de corriente (p$ = 
cte) y líneas equipotenciales o piezométricas ($ = cte). 


En estas condiciones tenemos: 


n/2a 
1 + (rrpp) (/a) cos(mp/2a) E 
A(x,y) — - |] (10) 
x/2a 
d1 + (órI-óI) (2/0) sen (19/20) 


(trabajamos en coordenadas polares (r, p) en lugar de cartesianas (x, y). 


Por otro lado elijamos la interpolación más simple posible, 
lineal: 


Y +11) 11 > Y 


2 e 2 ; 
Tr (€,n) = (1) 
1 + 11] A 11 - Y 
E 
2 2 


En consecuencia se deduce de (9) que la transformación 
geométrica necesaria, que reproduce exactamente la solución analítica (10) 
es: 


T=Altor (12) 


Desgraciadamente esta transformación no es cuadrática y por 
tanto no será posible reproducirla con elementos isoparamétricos 
cuadráticos. En la práctica esto significa que la aplicación T, en vez de 
aplicar, por ejemplo, la recta £ = 1 sobre la verdadera línea de corriente 
Y = Y, (línea de puntos de la figura 3) sólo la aproximará como se ob- 
serva en la misma figura, en el sentido de la interpolación, es decir hay 
coincidencia en los nodos (soporte). Con relación a la misma figura esto 
quiere decir que si nos fijamos en un punto P en el plano (£, y) tendremos 
que: 


n(P= (1, m-| $ ] 


Ahora bien, rm (P) es el valor interpolado en el punto P, ima- 
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gen de P por T, y queda claro que en P, y * py] porque la línea £ = 1 no 
coincide exactamente con la línea de corriente. Todo ello produce un 
error, que sin embargo se puede esperar que sea pequeño dado que en 
realidad verificamos (10) en un sentido medio. 


El resultado más interesante (y en cierto modo sorprendente) 
es que cuando a = Tm la expresión (12) se desarrolla sin excesiva 
dificultad a partir de (10) y (11), [5], para dar: 


b 
= (1+6)? - =(1:m)? 
4 4 x 
T(E,m) - - |] (13) 
ab y 
> (L+E) (den) 


Como quiera que la expresión (13) es cuadrática en (€,n) ele- 
mentos cuadráticos y de grado más elevados son capaces de reproducirla 
exactamente, y por tanto reproducir la solución analítica (10) alrededor 
de S sin ningún error. 


En el caso que nos ocupa, se trata de resolver un problema 
elástico y de aproximar por tanto los desplazamientos. basta para ello 
considerar el caso a = rm (y a=b sin pérdida de generalidad) y ubicar los 
nodos de acuerdo con la regla: 


x= + Ep? - (11? 
(14) 


E 2 
Yi = 2 Leg” 134 


Esto da lugar al elemento de la figura 4. 


Para obtener, por ejemplo, la posición del nodo 8, pondríamos 
en (14) £, >= -1, n, = 0, de donde x, = -a/4, y, = 0. Vemos que se trata por 
tanto de un elemento con nodos a 1/4 como los analizados en [3]. Ahora 


bien, este elemento tiene frente a aquéllos la ventaja de que 


(1) tiene lados curvos, lo que es importante para singularidades dis- 
tintas de r*, que aparecen en problemas elásticos cuando los 
dominios no son convexos, 


(11) es un elemento cuadrilateral, que se comporta de acuerdo a las 
fórmulas (10) que pueden aproximar convenientemente las (1) con una 
adecuada elección de los parámetros; además su comportamiento 
radial es el deseado como ya se ha demostrado en [18], y 
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(iii) permite modelizar una fractura con sólo 2 elementos, en contra de 
los 4 a 8 recomendados en [2] para los elementos a 1/4 clásicos. 


4. UN_TEST NUMERICO 
En la figura 5 aparece el problema tratado. 


Se trata de una pieza en deformación plana sometida a una 
tracción uniforme en los extremos y con una grieta inclinada en el borde. 


Para la resolución del problema se efectuó en primer lugar un 
modelo con un mallado bastante grosero (12 elementos y 49 nodos) que 
aparece en la figura 6. 


Numéricamente el primer problema que se plantea es el mal 
condicionamiento de la matriz de rigidez del sistema. Esta cuestión es 
clásica en problemas de elasticidad y da lugar a que el modelo sea estruc- 
turalmente inestable. Para resolverlo se ha usado un método de resolución 
iterativo (gradiente conjugado) con precondicionamiento diagonal [19]. El 
problema en tales condiciones converge en 104 iteraciones. 


Tabla 1 : Resultados del modelo Fig. 6 


Teórico| Error Teórico| Error 


[20] (3) [20] (%) 











Por otro lado se ha usado un orden de integración 2 x 2 para 
los cuadriláteros y 3 puntos para los triángulos, como es recomendable 
para este tipo de elementos [2]. Los resultados obtenidos se recogen en la 
tabla 1, en ella aparecen los factores de forma K,/0./na y K,,/0/rma cal- 
culados mediante (4) en los cuatro nodos más próximos al vértice de 
grieta. Como se ve el error es apreciable. Aunque se ha empleado una 
relación 1/a = 0.04 (donde 1 es la longitud del elemento singular), que 
entre dentro de los valores recomendados por diversos autores, [1], [2], 
[9], para otros elementos singulares y de los límites de aplicabilidad de 
(4), se comprueba que se obtienen mejores valores para los nodos más 
alejados de S, 


Este erróneo comportamiento del modelo es debido a que la 
relación entre el tamaño del elemento singular y el de los elementos 
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adyacentes es muy pequeña, lo que origina una mala representación de los 
campos asintóticos (1) fuera del elemento singular. 


Al objeto de mejorar los resultados del modelo anterior, se 
han desarrollado sendos modelos más refinados en la zona singular, que 
aparecen en la Figura 7 y la Figura 8. 


En el primero se han utilizado los elementos descritos en el 
$3 del presente trabajo, mientras que en segundo los seis elementos que 
rodean el vértice de la singularidad son triángulos standard, triángulos 
1/4 [16], y elementos de tipo Akin [11]. 


Por otro lado se han empleado dos condiciones de contorno dis- 
tintas (Figura 9) que como se observa en los resultados son equivalentes. 


Tabla II: Resultados para los modelos de las figuras 7 y 8 

















Error (3%) 
E : E 
1.09 1.09 0.42 0.42 

1 ' 8.91 









Modelo 1 
(Standard) 







ET 26.40| 26.41 
1.33 1.31 0.70 0.75 
F11 23.46| 31.02 
Fi 4.14 4.12 
Modelo 3 
E 1.15 1 E5 0.56 0.56 
0.59 0.59 


Modelo 2 
(Akin) 










Modelo 4 
83 
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Tabla 111: Resultados K* para los modelos de las figuras 7 y 8 





Error (3%) 











Modelo 2 





(Fórmula 7) 





Fi 18.30 |18.30 
Modelo 3 


(6a,b) 


.98 





.20 l E l 
Fi] 11.43| 12.07 


Se puede comprobar, Tabla 11, que los resultados son superio- 
res en cálculo en modo I que en modo II. Como es bien sabido el modo de 
fractura predominante es el Il. De hecho, la fisura siempre tiende a 
propagarse en modo I, [21] y precisamente el hecho de que estén presentes 
ambos modos de carga al comienzo de la propagación hace que la fisura se 
desvíe (buscando el modo I). Precisamente por este motivo, las expresiones 
para los desplazamientos verticales (u, y por tanto K,) dan mucho mejores 
resultados que las de los horizontales (u, y por tanto K,,), dado que son 
de mayor magnitud y mucho más próximos a las expresiones asintóticas (1). 
De hecho un aumento del ángulo [Pf hasta alcanzar el modo 1 puro (Pf = 1/2) 








daría lugar a mucho mejores aproximaciones, como ha sido abundantemente 
testeado en la literatura [15]. 


La fórmula (6a) (K*,) mejora los resultados del modelo con 
elementos curvos ($3), mientras que no lo hace para el modelo con elemen- 
tos 1/4 (Tabla 111) esto es bido al hecho de que la relación 1/a es mucho 
más grande para este tipo de elementos y de forma evidente es más exacta 
la fórmula (4a) como también se demuestra numéricamente en [15]. No sucede 
lo mismo con las (6b) lo que indica claramente que para el cálculo de K,, 
es recomendable utilizar los desplazamientos horizontales más próximos al 
vértice de la grieta y un ratio 1/a algo más elevado que para el cálculo 
de K;. 


Observemos también que los elementos de tipo Akin, dan unos 
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resultados bastante deficientes, esto es debido fundamentalmente al hecho 
de que como se demuestra en [4] la representación de las deformaciones 
aproximadas por este tipo de elementos no es la correcta (ver también 
[17]). Mejoras importantes se pueden obtener para K,, en cualquier caso, 
con el empleo de (7). 


Digamos por fin que el uso de elementos normales conduce a 
errores elevados, como es sabido, en el campo de la mecánica de la frac- 
tura. 


Digamos por fín que todos los cálculos se han llevado a cabo 
en doble precisión sobre un VAX/8300, usando como código de cálculo el 
programa MINIMEF [22]. 


CONCLUSIONES 


(1) Para problemas de fractura en modo mixto parece muy interesante el 
uso de elementos finitos singulares, dado que con ellos se pueden 
calcular los factores ce intensificación de tensiones d..rectamen!.e 
sin necesidad de un postproceso complejo. 


(11) Dentro de los elementos finitos singulares son recomendables los 
elementos a 1/4 (triángulos o cuadrángulos de lados curvos del $3) 
por su sencillez y buenos resultados. 


(iii) Los problemas en modo mixto exigen mallados más refinados alrededor 
de la singularidad (no tan simples como para fracturas, en modo 1 
puro) e incluso el uso de métodos ¡iterativos precondicionados de 
resolución. 


(iv) Para los cálculos es necesario trabajar con las componentes de los 
desplazamientos en la dirección de la grieta. Los desplazamientos 
verticales, de los que se deduce K,, dan mejores resultados que los 
horizontales, de los que se deduce K,,. este hecho ya había sido 
observado [10] y [15] y es la razón de que en general los elementos 
singulares deban comprobarse en modo mixto (el más severo 
"numéricamente"). Obsérvese sin embargo que el modo Il, es el 
predominante y el mejor calculado. 


(v) La relación l1/a es importante a la hora de elegir el método de 
cálculo de K, y K,,. En general conviene una relación 1/a pequeña 
(<0.1) y usar las fórmulas (4) que tienen un significado físico más 
evidente (dado que calculan K, y K,, a partir del "salto" de los 
desplazamientos) [2] y [15]. 


(vi) Las condiciones de contorno no influyen siempre que estén bien 
modelizadas en el cálculo de los factores de intensificación. 


(vii) Como conclusión general podemos decir que el cálculo de K, y K,, en 


modo mixto es importante, dado que supone una etapa previa en 
cualquier criterio de fractura adoptado, cálculo del ángulo de 
propagación y estudio de estabilidad; por ello y dado que es 
difícil emitir recomendaciones generales sobre el uso de elementos 
singulares, parece necesario efectuar varios tests numéricos para 
cada elemento y un estudio detallado del modelo de elementos 
finitos a usar en cada caso concreto. 
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Fig. 1 - Ejes locales en el Fig. 2 - Cálculo de los F.1.T. 
vértice de grieta. 








Fig. 3 - Elemento singular 








Fig. 4 - Caso a = “um 
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Fig. 5 - Problema estudiado Fig. 6 - Modelo grosero 








Fig. 7 Modelo 4 Fig. 8 - Modelos 1, 2 y 3 
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Fig. 9 - Condiciones de contorno 
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